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Abstract 


In this study, putative interactions between all of the retinoic acid (RA) ligands (all-trans (At), 9-cis (9c), 
and 13-cis (13c)), and VEGF receptors (VEGFR-1, -2 and -3) were investigated. It was performed considering the 
glycosylation status of the receptors to achieve a more reliable mode of interactions based on glycomics. We found 
that RAs may have a higher affinity for ligand-binding domains in VEGFRs. Furthermore, all RA isomers can 
strongly attach to VEGFR-3 receptor in comparison to other ones. It was also demonstrated that receptor dimerization 
of RAs may be less targeted. Moreover, regarding post-translational modifications, glycosylated structures showed 
conflicting binding energies. RAs may target the human vasculature, specifically lymph vessels, through VEGFR-3. 
In addition, the ligand binding-mediated activation of VEGFRs may be affected by these agents. Also, the 
glycosylation status of the receptors can interfere with these manners. Furthermore, our results confirmed that the 
consideration of carbohydrates in crystal structures is essential for a better interpretation of ligand/receptor 
interactions during drug discovery studies. Even though these observations improved our understanding of the 
binding patterns of RAs to VEGFRs, validation of these results needs further analysis to introduce these biomolecules 


as anti-VEGF remedies. 


Keywords: Retinoic acids, Angiogenesis, VEGF receptors, Glycomics, in silico 


Introduction 


Organization of blood and_ lymphatic 
vessels, vasculogenesis and angiogenesis, is critical 
to provide whole-body with fresh oxygen and 
nutrient supply and remove catabolites. At the 
cellular level, activation of vascular endothelial 
growth factor receptors (VEGFR-1, -2, and -3) by 
their cognate ligands (different isotypes of VEGF) is 
the most important receptor tyrosine kinase/RTK- 
related signaling pathway. This pathway adjusts 
multiple processes which are essential for 
developmental, physiological, and pathological 
neovascularization (Christensen et al. 2017; Lee et 
al. 2017; Qiao et al. 2006; Shibuya 2013). The 
VEGFs are signal proteins which are produced by 


elimination of both pro-peptides, they obtain binding 
affinity for WEGFR-2. VEGF-E which exists in 
poxviruses can specifically bind to WEGFR-2 
(Mercer et al. 2002). 

As the ligand binding to the extracellular domain 
(ECD) of VEGFRs is required for their stimulation, 
the disruption of these ligand/receptor complexes 
deregulates their RTK and normal physiological 
activities. For this reason, ligand trapping from 
VEGFR-1 is resulted in the placental loss of 
preeclampsia and abnormalities in retinal and 
corneal vascularization (Markovic-Mueller et al. 
2017). In contrast, VEGFR-1 traps the main 
angiogenic ligand VEGF-A from VEGFR-2 and 
displays a negative regulatory role against VEGFR- 
2-induced blood vessel sprouting. Alternatively, 


cells to trigger the formation of vessel networks via VEGFR-2/VEGFR-3 heterodimerization is an 
binding to their corresponding RTKs. All isoforms additional signal to facilitate this sprouting. 
of VEGF-A can bind to VEGFR-1 and VEGFR-2, Moreover, VEGFR-3 homodimers mediate 


while VEGF-B is special for VEGFR-1 (Alvarez- 
Aznar et al. 2017; Goel and Mercurio 2013; Jeltsch 
et al. 2006; Leppanen et al. 2010). VEGF-C and 
VEGE-D isoforms with cleaved C-terminal domain 
are high-affinity ligands for WVEGFR-3. Upon 
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lymphogenesis. So, a defective VEGFR-3 cascade 
results embryonic death (Leppdnen et al. 2013; Qiao 
et al. 2006). 

Considering the fact that classical anti-VEGF 
compounds are clinically administrated to prevent 


Licensed as Open Access CC-BY 


Journal of Cell and Molecular Research (2021) 12 (2), 87-96 


particular disease (cancer, retinopathy, 
endometriosis, and Xe) forth)-associated 
neovascularization, the discovery of detailed 


molecular mechanisms of blood vessel formation is 
essential to design advanced anti-angiogenic 
remedies (Hegde et al. 2018; Ma & Ning 2019). 

The glycome is the entire complement of sugars, 
whether free or present in more complex molecules 
of an organism. Glycomics is a comprehensive and 
multi-aspect study of glycome (Rudd et al. 2017). 32 
types of sugar linkages were reported in various 
Saccharides. These building blocks can increase the 
degree of complexity. Sugar structures, are highly 
branched. Glycans are highly dynamic as they can 
bound to proteins or conjugate with lipids to form 
modified structures (Aizpurua-Olaizola et al. 2018). 
Glycans play significant roles during viral/bacterial 
recognition, cell signaling events, intrinsic immunity 
modulation, prohibition of cell proliferation, cancer 
expansion, cell fate determination, invasion, 
circulation arrangement, protein folding and other 
important biological procedures (Aizpurua-Olaizola 
et al. 2018; Rudd et al. 2017). Hence, glycomics is 
widely applicable in the clinic as it was reported 
previously (Reid et al. 2012). 

Currently, we know that glycosylation facilitates 
signal transfer from the surrounding 
microenvironment to inside the cell. The structure of 
carbohydrates has direct impact on ligand/receptor 
binding, dimerization, and internalization of various 
signaling complexes, and, thus activation of 
different cell receptors. RTKs are the most prevalent 
type of Asparagine (N)-linked glycans that transduce 
signals which are important for cell 
communications. Accordingly, it could be 
concluded that VEGFRs mediate main functional 
properties of the vascular endothelial cells during 
angiogenesis. There are 13, 18, and 12 N- 
glycosylation sites with different carbohydrate side- 
chains in VEGFR-1, -2, and -3, respectively 
(Contessa et al. 2010; Itkonen and Mills 2013; 
Lopez-Sambrooks et al. 2016). Although it has not 
been fully understood so far, these sugar residues, 
which are not completely detected in crystal 
structures, may also affect the conformational 
features of VEGFRs. Therefore, comprehending the 
effects of the sugar structures on the binding of a 
candidate drug VEGFRs improves our knowledge 
about the binding mode of action of drug/VEGFRs 
such complexes. also, this direct impact is important 
for a glycomics-based point of view in drug design 
against VEGFRs. 


Retinoids are physiological derivatives of 
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vitamin A which have several positive impacts on 
normal cell growth and differentiation, tissue 
homeostasis, developmental organogenesis, and 
visual performance (Khalil et al. 2017; Mallipattu 
and He 2015; Zhu et al. 2015). However, they may 
have destructive effects on the morphogenesis 
during human embryogenesis (Gudas, 1994; R6nn et 
al., 2015). All retinoids including all-trans (At), 9- 
cis (9c), and 13-cis (13c) retinoic acids (RAs) 
emerge more recently as antineoplastic agents. 
Accordingly, on the contrary to the harmful 
influence of At-RA, a main active isomer of RA, 
(Rihl et al. 2018; Tsuji et al. 2015; Hu et al., 2020) 
that frequently causes embryonic death due to 
developmental deformity, it has been used in cancer 
therapy towards acute promyelocytic leukemia 
(Huang et al., 1988) colorectal cancer, pancreatic 
cancer, oral leukoplakia and skin cancer (Applegate 
and Lane 2015; Di Masi et al. 2015; Lodi et al. 2016; 
Moon et al. 1997; Tarapcsak et al. 2017; Uray et al. 
2016; Szymanski et al., 2020). 

Interestingly, although we know that RAs are 

diffusible biomolecules (Minkina et al. 2017) and 
exert their impressions on gene regulation via 
nuclear receptor RAR/RXR. However, several 
studies demonstrated that RAs can tune-up the 
angiogenesis through binding to VEGFRs with an 
unclear pattern (Njar et al. 2019; Simandi et al. 2018; 
Urushitani et al. 2018; Costantini et al., 2020). 
In the present study, via a computational-based 
approach we investigated the binding manner of At- 
RA, 9c-RA, and 13c-RA with the extracellular 
regions of three human VEGFRs (R-1, R-2, and R- 
3). We found an alternative mechanism which is 
responsible for the VEGFR-mediated inhibitory 
implications of RAs_ during angiogenesis. 
Furthermore, regarding the importance of VEGFRs’ 
saccharide residues inactivation, we also compared 
the possible binding sites of three isoforms of RAs 
to both glycosylated and unglycosylated crystal 
structures of VEGFRs introduced so far. 


Materials and Methods 


Preparation of structures 

Retinoic acids are small biomolecules with an 
average mass of about 300.435 kDa. The structures 
of interest isoforms (At-RA, 9c-RA, and 13c-RA) 
have been retrieved from PubChem database 
(http://pubchem.ncbi.nlm.nih.gov/) in structure- 
Data file (SDF) format which is an acceptable input 
for our docking program. Identifier (ID)s_ of 
retrieved isoforms from public repository were as 
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follows: 6419707; 449171, and 5282379, 
respectively (Dunning 2018). Also, the PDB files of 
human VEGFR1-3 have been obtained from the 
Protein Database  (http://www.pdb.org). — All 
structures were prepared _ by WebLab 
(http://weblab.cbi.pku.edu.cn) (Liu et al., 2009) and 
the visual molecular dynamics (VMD 1.9.3) package 
(www.ks.uiuc.edu/Research/vmd/) (Zhang 2015). 
The Swiss-PdbViewer modeling program was 
applied to perform energy minimization of all 
investigated structures 
(http://www.expasy.org/swissmod) (Guex et al. 
2009). As the cut-off for bond, torsion, angle, 
improper angles, non-bounded and _ electrostatic 
functions was set for 10 A by default, energy (E) 
exchanges between two steps was stopped below 
0.05 kJ/mol, and the acting forces between any atom 
stopped below 10 A. 
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The structural details including PDB IDs, domains, 
length, and glycosylation status of the extracellular 
regions of our studied VEGF receptors are 
summarized in the Table 1. 

In order to minimize the energy of all 3D structures 
prepared by the Web Lab, we performed the energy 
minimization procedure using Swiss-PdbViewer 
(DeepView). The first (E1) and the end (E2) energy 
amounts for each molecule were listed in the Table 
2. Moreover, docking energies obtained for different 
RA-VEGFR complexes using HEX 8.0 were 
described in Table 1. 

A comparison of assessed docking energies for 
various RA-VEGFR complexes in this work 
indicated that At-RA, 9c-RA, and 13c-RA (-256.10, 
-256.11, and -248.56 kJ/mol, respectively) are 
highly bound to VEGFR-3 4BSK (Table 1). This is 
also occurred at different amino acid residues in 
VEGFR-3 D2 (OA-133: ARG-CA to OA-190: 


Conaeateentleanae! 


A) At-Retinoic acid (At-RA) 


B) 9-cis-Retinoic acid (9c-RA) 


COOH 


C) 13-cis-Retinoic acid (13c-RA) 


Figure 1. Schematic structural representation of retinoic acid isoforms including all-trans (At), 9-cis (9c), and 13- 
cis (13c) retinoic acids (RAs). At-RA is the main active isomer of RA (http://pubchem.ncbi.nlm.nih.gov/). 


Docking parameters 

Docking studies were carried out using HEX 8.0 
(http://hexserver.loria.fr), which is an interactive 
molecular graphics program to calculate and display 
feasible docking modes of pairs of protein and ligand 
molecules (Roudini et al. 2020). Moreover, it can 
predict the position of ligand-receptor contacts 
approximately. All docking parameters were left as 
their default values (Pourhashem et al. 2017). 


Results 


The molecular structures of retinoic acid isoforms 
including At-RA, 9c-RA, and 13c-RA_ were 
represented in Figure 1 as they were obtained from 
PubChem. Chemical differences and reactive groups 
are shown for each structure. 
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GLY-CA) (Figure 2). 

Besides, our ligands can attach to other 
investigated structures with lower tendencies and 
they have the lowest binding mode in complex with 
VEGFR-1 D2 (1QTY), as demonstrated in Table 1. 
Also, we found At-RA, 9c-RA, and 13c-RA bind 
slightly to VEGFR-2 and VEGFR -1 D3 (residues 
from [O] R-143: ASN-CA to [O] R-278: LYS-CA 
and [O] X-241: LYS-CA to [O] X-331: LYS-CA, 
respectively) (Figure 2). In addition, as it was 
demonstrated in Figure 3 regarding the trends of 
binding energies for RA-VEGFR complexes, a 
decreased binding mode was detected for VEGFR-1 
D3 (ST89) in comparison to VEGFR-2 D3 (2X1X) 
and also for VEGFR-2 D3 (2X1X) in comparison to 
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able 1. Structural details ain, Domain, Glycosylation status) for s and docking simulations for S- 
VEGEFRs. 
Receptors PDB Chain Domain Glycosylation status Docking Energy (kJ/mol) 
IDs (amino acid Number of At-RA 9c-RA  13c-RA 
positions) Sugars and aa positions glycosylation 
sites 
VEGFR1 1FLT 95aa D2 - - -228.38  -223.9  -209.82 
(Phe X135- (aal51— 
Thr X226) aa 214) 
1QTY 10laa D2 - - -15.76 -17.29 -18.19 
(Phe T135- (aal51 — 
Gln T225) aa 214) 
S5T89 646aa D1 ASN X 100 (NAG X 701 -217.79 -231.50  -218.01 
(Asp 31-Thr (aa32-aa123), & NAG X 702), ASN X 
Y654) D2 196 (NAG X 703), ASN X 
(aal51-aa214), 251 (NAG X 704), ASN X 
D3 323 (NAG X 705), ASN X 
(aa230-aa327), 402 (NAG X 706), ASN X 
> 417 (NAG X 707 & NAG 
D4 X 710), 
(aa335-aa421), Asn x 574 (NAG X 708), 16 
D5 ASN X 625 (NAG X 709) 
(aa428-aa553), / ASN Y 100 (NAG Y 
D6 701), ASN Y 164 (NAG Y 
(aa556-aa 654) 702), ASN Y 196 (NAG Y 
703 & NAG Y 704), ASN 
Y 251 (NAG Y 705), ASN 
Y 402 (NAG Y 707), ASN 
Y 417 (NAG Y 708), ASN 
Y 547 (NAG Y 709), ASN 
Y 625(NAG Y 710) 
VEGFR2 2X1W 213aa D2 ASN L 143 (NAG L -244.06 -244.45 -241.28 
(Phe L125- (aa 41-aa207), 2401), ASN L 245 (NAG 
Glu L 326) D3 L 2601 & NAG L 2602), 
(aa224-2a320) ASN L 318 (NAG L 
2701)/ ASN M 245 (NAG 
M 1201 & NAGM 1201), 
ASN M 318 (NAG M 
1301, NAG M 1302 & 12 
BMA M 1303)/ ASN N 
143 (NAG N 1401), ASN 
N 245 (NAG N 1601), 
ASN N 318 (NAGN 1701 
& NAG N 1702)/ ASN O 
143 (NAG O 2001), ASN 
O 158 (NAG O 2101), 
ASN O 245 (NAG O 2201 
& NAG O 2202), ASN O 
318 (NAG O 2301). 
2X1X 213aa D2 ASN R 143 (NAG R 404), -251.99 -245.43  -242.62 
(ProR124-  — (aal41-aa207), ASN R 245 (NAG R 402 3 
Glu R326) D3 & NAG R 403), ASN R 
(aa224-aa320) 318 (NAG R 405 & NAG 
R 406). 
3KVQ 108aa D7 - : -241.31 -247.96  -240.43 
(Ile A669- (aa667-aa753) 
Gly A756) 
3V6B 424aa D2 - - -189.43  -206.67 = - 186.63 
(His R131- (aal41-aa207), 
Phe R329) D3 
(aa224-aa320) 
VEGFR3 4BSJ 232aa D4 ASN A 411 (NAG A 601), -51.50 -49.20 -43.38 
(Asp A328- — (aa331-aa415), © ASN A515 (NAG A 602) 2 
His A554) D5 
(aa422-aa552) 
4BSK 214aa D1 ASN A 33 (NAG A 301), -256.10 -256.11 -248.56 
(Pro A30-He — (aa30-aal27), | ASN A 104 (NAG A 302), 3 
A224) D2 ASN A 166 (NAG A 303) 


(aal51-aa213) 


*ASN: Asparagine 


* NAG: N-acetyl glucose 
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Table 2. Energy minimization (first and end energy) for VEGFRs. 


Receptor PDBIDs Chain Fi (first energy minimization E2 (end energy minimization 
kJ/mol) kJ/mol) 
VEGFRI1  1FLT X, Y -6040.213 -10363.84 
1QTY T, U, X, Y -4588.174 -22151.435 
ST89 XxX, Y -40754.605 -71025.164 
VEGFR2 2X1W L, M,N,O -23284.633 -37331.891 
2X1X R +343794.188 -11717.921 
3KVQ A -1688.371 - 4945.834 
3V6B R -2886.585 -10666.744 
VEGFR3 4BSK A 16982.123 -9005.400 
4BSJ A -7782.256 -12937.562 


Al) AtRA 


Psat 


: 


! 


Bl) AtRA 


x 


\ 


| 


Cl) At-RA 


C2) 9c-RA C€3)13¢-RA 


Figure 2. Docking of three retinoid types (At-RA, 9c-RA, and 13c-RA) with VEGFR-1, 2, and 3 complexes: A1-3, 
B1-3, and C1-3 are docking results of retinoid types with VEGFR-1 (PDB-ID 5T89), VEGFR-2 (PDB-ID 2X1X), and 
VEGFR-3 (PDB-ID 4BSK), respectively. RA isoforms differentially bind to different domains of VEGFRs. 
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5T89 | 2X1X | 4BSK | 


At-RA 


5789 | 2X1X | 4BSK | 


5T89 | 2X1X | 4BSK 


9c-RA 13-RA 


Figure 3. Docking energy amounts are shown for various complexes of At-RA, 9c-RA, and 13c-RA isoforms with 
different VEGF receptors, including VEGFR-1 D3 (5T89), VEGFR-2 D3 (2X1X), and VEGFR-2 D3 (2X1X). 


VEGFR-3 D2 (4BSK) in a gradual mode in case of 
9cRA isoform. A similar trend was observed for both 
At-RA and 13-RA isoforms, while differences were 
more remarkable for 5T89 in comparison to 2X1X 
and 4BSK. According to crystal structures of three 
VEGFRs, structures 5T89, 2X1W, and 4BSK, for 
VEGFR-1, -2, and -3, have resolved the higher 
number of carbohydrate residues among other 
structures studded up-to-now, respectively (Table 
1). Considering the glycosylated WEGFR-2 
structures, 2X1W (with a higher glycosylation 
status) has a lower docking tendency (At-RA: - 
244.06; 9c-RA: -244.45; 13c-RA: -241.28, kJ/mol) 
in comparison to a less glycosylated structure 2X1X 
(At-RA: -251.99; 9c-RA: -245.43; 13c-RA: -242.62, 
kJ/mol). These observations were not confirmed for 
other receptors such as glycosylated 2X1X with a 
lower energy versus non-glycosylated 3V6B with a 
higher energy (Table 1). 


Discussion 

The regulatory function of the VEGF signaling 
system has comprehensively been established in 
normal de novo and pathological revascularization. 
This cascade is modulated by the activation of 
VEGF receptor family members (VEGFR-1, -2, and 
-3). VEGFR-2 triggers the blood vessels 
organization whereas VEGFR-1 plays a negative 
feedback for VEGFR-2 function (Greenberg et al. 
2008; Stuttfeld and Ballmer-Hofer 2009). In 
VEGFR-3-related signaling pathways, VEGFR-3 
(only) induces lymphatic vascularization and also in 
co-operation with VEGFR-2 mediates blood 
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angiogenesis (Leppdnen et al. 2013). These trans- 
membrane RTKs_ have’ seven extra-cellular 
immunoglobulins (Ig)-like domains (D1-7) and their 
activation process involves the ligand binding 
through D1-3, dimerization through D4-7, and auto- 
phosphorylation of the cytoplasmic kinase domains. 
The ligand-mediated induction mode via a 
membrane-distal domain 2 (D2) is conserved for all 
VEGF receptors. Additionally, for both VEGFR-1 
and 2, D3 (not mainly) is used for ligand attachment 
(Leppanen et al. 2010; Leppanen et al. 2013; 
Markovic-Mueller et al. 2017). 

Recent studies disclose the potential of retinoids 
to interrupt with angiogenic processes. Weninger et 
al. in 1993 illustrated that retinoids are potent 
inhibitors of VEGF/VPF production by normal 
human keratinocytes. These cells may contribute to 
the therapeutic effects of retinoids (Weninger et al. 
1998). Also, in 2007 Noonan et al. showed that the 
synthetic retinoid 4-hydroxy fenretinide (4HPR) 
displays anti-angiogenic effects (Noonan et al. 
2007). In addition, microarray data by Albini et al. 
demonstrated that some anti-angiogenic agents, such 
as N-acetyl-cysteine and 9c-RA, induce molecular 
indices in endothelial cells which mimic in vitro 
senescence (Albini et al. 2012). Tsuji et al. showed 
evident data for the binding of retinoic acid isomers 
(At-RA, 9c-RA, and 13c-RA) with RXRs and RARs 
(Tsuji et al. 2015). Moreover, Zhong et al. confirmed 
that CYP26C1 is a 4-oxo-atRA hydroxylase and 
probably be important for adjusting the condensation 
of active retinoids in human tissues (Zhong et al. 
2018). 

Taken together, there is no clear data indicating 
that how RAs affect VEGFRs from outside the cell. 
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Here for the first time, the binding patterns of RAs 
with VEGERs were clarified using a docking-based 
approach. Results revealed that all three RA isomers 
(At-, 9c-, and 13c-RA) can remarkably bound to 
VEGFR-3. In comparison to interactions with other 
domains, RAs have a higher affinity to D2 in 
VEGFR-3 (Figure 2). This may suggest a novel 
mechanism for these agents which targets the 
vasculature, especially through the ligand (VEGF-C) 
detachment from its respective receptor which leads 
to the lymph/blood vessels deterioration. These 
biomolecules can also interfere with blood vessels 
formation through the inhibition of WEGF-A- 
VEGFR-2 and (with a lower affinity) VEGFR-1 
(Table 1). Although RAs may have strong contacts 
with similar domains, as depicted in VEGFR-3 D2 
(Figure 2) or VEGFR-1/-2 D3 (Figure 2), they can't 
interact with similar residues, may be owing to 
variations in RA isomers (Figure 1). Actually, 
structural remodeling for receptor dimerization in 
D4-7 could not be sharply affected by RAs, 
considering the competed energies which were 
demonstrated in Table 1. As the kinase activity of 
VEGFRs is fundamentally concerted by their 
glycosylation status (Gomes Ferreira et al. 2018), 
their glycosylated 3D_ structures were also 
considered in the current study. In this regard, 
glycosylation may lead to a lower affinity of RAs for 
VEGFR-2 (2X1W) when compared with 2X1X. 
Conversely, glycosylation may affect the ST89 
folding in a way that different modes achieve in RA- 
VEGFR-1 docking assessments. These 
inconsistencies in binding energies between 
glycosylated and unglycosylated structures suggest 
that sugar residues can alter the ligand-binding 
manners. Accordingly, post-translational 
modifications like glycosylation should be 
considered in such docking studies to reach a reliable 
behavior. It implies the importance of glycomics- 
based approaches in ligand/receptor binding studies. 

Taking a computational analysis, instead of a real 
experiment, is easy and cost-benefit, though, 
validation of the outputs by real experiments are 
inevitable. It is clear that all factors applied during 
simulations are based on approximation and are 
under control by users. Currently, by advances in 
computational biology docking of several thousands 
of ligands were performed which can be applied in 
pharmaceutical industries to design improved 
structure-dependent drugs (Yadav et al. 2018). 
Nevertheless, we can have a better interpretation of 
the applications of RAs in tissue regeneration, 
wound repair, and cancer therapy or other VEGF- 


http://jcmr.um.ac.ir 


Fallah et al. 


associated system uses in the future. thus, molecular 
dynamics and experimental analysis in further 
investigations may confirm our findings. 
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